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We present the time-dependent restricted-active-space self-consistent field (TD-RASSCF) theory 
as a new framework for the time-dependent many-electron problem. The theory generalizes the mul- 
ticonfigurational time-dependent Hartree-Fock (MCTDHF) theory by incorporating the restricted- 
active-space scheme well known in time-independent quantum chemistry. Optimization of the or- 
bitals as well as the expansion coefficients at each time step makes it possible to construct the 
wave function accurately while using only a relatively small number of electronic configurations. 
In numerical calculations of high-order harmonic generation spectra of a one-dimensional model 
of atomic beryllium interacting with a strong laser pulse, the TD-RASSCF method is reasonably 
accurate while largely reducing the computational complexity. The TD-RASSCF method has the 
potential to treat large atoms and molecules beyond the capability of the MCTDHF method. 

PACS numbers: 31.15.-p, 32.30.-r, 33.20.Xx 
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I. INTRODUCTION 

Development of reliable approximate theories for the 
description of time-dependent many-electron dynamics 
has been desirable for decades, and its importance is es- 
pecially emphasized nowadays by the need of support 
to experiments on the real time analysis and control 
of ultrafast electronic and nuclear dynamics of atoms 
and molecules by intense laser pulses [1-5]. In numer- 
ical simulations, however, the problems are most often 
simplified by the single-active-electron (SAE) approxi- 
mation [6] , which assumes that only one electron is mov- 
ing in an effective potential. In theoretical approaches, 
the combination of the SAE and the strong field ap- 
proximations, where the interaction with the atomic or 
molecular potential is treated perturbatively, has been 
widely accepted as a standard approach in this research 
area. The Lewenstein model [7], which is built on these 
assumptions, makes it possible to easily compute high- 
order harmonic generation (HHG) spectra of atoms and 
molecules and also provides a clear physical picture based 
on the semiclassical three-step model [8-10]. While 
the studies within the framework of the SAE approx- 
imation have succeeded in providing a simple picture 
of the phenomena, multielectron effects are also recog- 
nized to play a crucial role, e.g., in time delay studies 
of photoionization [11, 12], and moreover multiple or- 
bital contributions to HHG spectra are widely observed 
for atoms and molecules [13-16]. To describe many- 
electron dynamics, several ah initio approaches beyond 
the SAE approximation have been developed. Among 
others, the time-dependent i?-matrix method is one of 
the most elaborate ways for describing single electron 
ionization processes and taking into account the elec- 
tron correlation [17-19]. One of the computationally and 
conceptually simpler approaches is the time-dependent 
configuration-interaction singles (TD-CIS) method [20- 
24] , in which the CTexpansion is truncated at singly ex- 
cited configurations relative to the Hartree-Fock ground 



state. Both these methods can be considered to be 
special cases of a more generalized concept, namely, 
the time-dependent restricted-active-space configuration- 
interaction (TD-RASCI) method [25]. 

Over the last decade, originating from the time- 
dependent Hartree-Fock theory [26-28], a more so- 
phisticate framework called multiconfigurational time- 
dependent Hartree-Fock (MCTDHF) theory has been de- 
veloped and quite recently shown its potential for ana- 
lyzing ultrafast laser driven electron dynamics in atoms 
and molecules [14-16, 29-35], and moreover for elucidat- 
ing the role of electron-nuclear correlation in a molecule 
during ionization [36-40] (see also references on the mul- 
ticonfigurational time-dependent Hartree (MCTDH) the- 
ory, e.g., the original paper [41], a review [42], and a 
textbook [43]). In the MCTDHF theory with the spin 
restricted ansatz, the A'o-electron wave function is ex- 
pressed by 



\^{t))= J2 ciimAt)), 

/GVfci 
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where |"I>/(i)) denotes a normalized Slater determinant 
for Nc electrons built from a set of time-dependent ac- 
tive spin-orbitals {4>i{t)}'i^i, and the multi-index / = 
(«i, • • • ,«jVc) is a string of orbital indices of which the 
Fock space is composed: Vpci = {(*!;••• ;*Afe)|l ^ «i < 
••• < iNa ^ 2M}. Using the Dirac-Frenkel-McLachlan 
time-dependent variational principle [44-47], the equa- 
tions of motion are derived which optimize the orbitals 
as well as the expansion coefficients in each time step. 
This optimization procedure leads to the expectation 
that the system can be accurately described with a rel- 
atively small number of orbitals, 2AI. However, because 
the Fock space Vfci in the MCTDHF theory is spanned 
by all possible configurations for a given set of spin- 
orbitals, the computational cost due to the expansion 
coefficients {C/j/gVpci is proportional to the number of 
ways that iVo electrons can be distributed in the 2M spin- 



orbitals 
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0((2M) 



Na 
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i.e., roughly speaking, exponentially scaling with respect 
to the number of electrons, N^. Hence for the investiga- 
tion of the nonperturbative laser driven electron dynam- 
ics, this unfavorable scaling with Nc impedes the MCT- 
DHF method to be extended to systems with more than 
a few electrons, i.e., beyond systems like He [30], Be [31], 
H2 [32, 33, 36], and LiH [34, 35, 38]. In order to carry on 
the study for larger systems, it is therefore inevitable to 
abandon the full Cl-expansion. 

For the time-dependent problem, it is natural to follow 
standard approximations in quantum chemistry such as, 
for example, CIS and explore their time-dependent coun- 
terparts, the TD-CIS method [20-24]. It is thus natural 
as well to investigate the possibility of a truncation of 
the expansion in the MCTDHF theory at a specific ex- 
citation level. Within the framework of the MCTDH 
theory, the truncation has already been explored [48]. 
Our study is hence aiming at a further generalization of 
the MCTDHF theory by incorporating it with the RAS 
approach; decomposing the single-particle Hilbert space 
into several subspaces, among which electron transitions 
are allowed with several restrictions. The new framework 
is hereafter referred to as the time-dependent restricted- 
active-space self-consistent field (TD-RASSCF) theory. 
It is emphasized that key ingredients of the theory are 
(i) use of time-dependent orbitals and (ii) truncation of 
the Cl-expansion. Despite of its simple concept, the trun- 
cation of the expansion partly destroys the principal bun- 
dle structure inherent in the MCTDHF theory [47, 49]. 
Accomplishing the truncation hence requires a careful 
analysis of the structure of the theory. A consistent for- 
mulation of the TD-RASSCF theory is the main purpose 
of the present work. 

This study is inspired by a recent work [50] , where an- 
other new method called orbital adaptive time-dependent 
coupled-cluster (OATDCC) theory was formulated. In 
this method, the CC-expansion is truncated at doubly 
excited configurations, but higher excitations are also 
partly included due to the nonlinear property inherent 
in the CC ansatz (see, e.g., Refs. [51, 52] for a discussion 
of time-independent CC theory). Furthermore, the CC- 
expansion ensures the size-consistency and -extensivity, 
which will be of utmost importance for correctly de- 
scribing dissociation processes of molecules. However, 
there are still some problems remaining: because the 
left- and right-wave functions in the OATDCC theory 
are parametrized in different ways, imaginary time re- 
laxation is not readily feasible as a means to calculate 
the ground state wave function needed for the following 
real time analysis of the dynamics. It is thus attractive 
to develop methods not suffering from these complica- 
tions while only slightly compromising the accuracy; the 
TD-RASSCF theory is one such example. 

The paper is organized as follows. The TD-RASSCF 
theory is formulated in Sec. II. The working equations 



are derived based on the time-dependent variational prin- 
ciple combined with the Lagrange multiplier method. 
Explicit forms of the equations are given and compared to 
the corresponding ones in the MCTDHF theory. A cen- 
tral aspect of the formulation is given in Sec. HI, where 
we concentrate on a discussion of the parametric redun- 
dancy in the time-dependent SCF theory. As a proof-of- 
principal application of the theory, one-dimensional (ID) 
model atoms are investigated in Sec. IV: calculations of 
the ground state wave function, followed by computations 
of the HHG spectra. The analysis of the convergence 
property of the TD-RASSCF calculations also uncovers 
the time-dependent many-electron dynamics. Section V 
provides a summary and concludes. A discussion of or- 
bital rotations and the accompanying simplifications in 
the case of a two-electron system is deferred to Appendix 
A. 



II. 



FORMULATION 



Consider an A'o-electron system described by a generic 
time-dependent Hamiltonian including one- and two- 
body operators. In second quantization, the Hamiltonian 
reads 



H{t) 



E 



KiU 



p^q 



1 



Vqlit)cWrCsCq, (3) 



where Cp (cD is the annihilation (creation) operator of 
an electron in the spin-orbital \(f>p{t)), satisfying the an- 
ticommutation relation {cp, cj} = 6^. The prefactors of 
the operators in the Hamiltonian read in first quantiza- 
tion: 



Kit) 



4il{z, t)h{r, t)(j)g{z, t)dz, 



(4) 



and 



v^/At) ^ J 4i^i,t)4{^^2,t) 

xv{ri,r2)(f)q{zi,t)(f)s{z2,t)dzidz2, (5) 

where the spin-orbitals are represented in the spin and 
spatial coordinates z = {r,a). The two-body opera- 
tor v{ri,r2) denotes the Coulomb repulsion between two 
electrons, and the one-body operator h{r, t) depends on 
time via dipole interactions with external laser fields. In 
the TD-RASSCF theory, the A'c-electron wave function 
is expressed by 



i*w>= E ci{mi{t)) 



(6) 
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in which, differently from Eq. (1), the Fock space Vras is 
now composed of several selected configurations and not 
the full configuration space. In this study, the orbitals are 
supposed to be classified as in Fig. 1: the single-particle 
Hilbert space is divided into two subspaces, T'-space con- 
tributing to the construction of the wave function and 



the supplementary virtual orbital space hereafter called 
Q-space. The indices p,q,r,s ■ ■ ■ denote orbitals belong- 
ing to either space, while the T'-space orbitals are la- 
beled by i, J, fc, Z, • • • , and the virtual Q-space orbitals by 
a, 6, c, d, • • • . The RAS scheme is based on dividing the 
T'-space into three subspaces usually denoted by RASl, 
2, and 3. In the original and the most general definition, 
the RASl and 3 spaces are characterized by the minimal 
and maximal occupation numbers, respectively, and the 
RAS2 space has no constraint [52-54] . The RAS scheme 
in this paper is, however, supposing more specific cases 
as shown in Fig. 1. Here the P-space is composed of an 
inactive-core space, Vq, and two active spaces, Vi and 
7^2, between which electron transitions are allowed with 
several restrictions (see Sec. IV). 

For describing the dynamics within this framework, we 
need the set of equations obeyed by the expansion coef- 
ficients and the orbitals. To this end, we follow a stan- 
dard prescription and use the Dirac-Frenkcl-McLachlan 
time-dependent variational principle [44-47] . Henceforth 
the explicit time-dependence of the parameters and the 
operators is dropped for brevity as long as there is no 
confusion. First we define an action functional (atomic 
units are used throughout) 

S[{CJ},{cf>^},{e)}]^ 




SI 



dt. 



(7) 
Then we use that a stationary point 

6S[{Cj},{<t>^},{e)}]=0 (8) 

provides the best approximation of the dynamics for the 
given ansatz. Here £^- is the Lagrange multiplier that 
ensures orthonormality of the "P-space orbitals during the 
time interval [0, T]. The variation of the action functional 
is taken as follows: 



SS[{CJ},{cf>^},{e)}] 



{S^\lt--H\\^) + {^\ 
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FIG. 1. Illustration of the division of the single-particle 
Hilbert space in the TD-RASSCF theory. The wave func- 
tion is expand by using orbitals in P-space, which is com- 
posed by an inactive-core space, Vo, and two active spaces, 
Vi and P2, where a partition exists through which electrons 
can transit with several restrictions. In accordance with the 
convention in the MCTDHF theory, the rest of the single- 
particle Hilbert space spanned by the virtual orbitals is re- 
ferred to as Q-space. The orbitals in either V- or Q-space are 
labeled by p,q,r,s,- ■ ■ , while the P-space orbitals are labeled 
hj i,j,k,l,- ■ ■ , and the Q-space orbitals by a,b,c,d,- ■ ■ . 



which ensures the orthonormality of the P-space orbitals 
at all times. Stationary conditions with respect to small 
variations of the other parameters 



65/ SC^ =SS/{S(f),\ =0 



(12) 



give us the equations of motion. Since the left- and right- 
wave functions are hermitian conjugates of each other, 
the stationary conditions SS/SCi = SS/\6(l>i) = result 
in a set of equations which is the hermitian conjugate of 
the set obtained from Eq. (12). 



E^4 



[{<i^t\4>j 



dt, 



(9) 



where, performing integration by parts, a time-derivative 
operator with a leftward-arrow is introduced to denote 
its action on the bra- vector. The variation of the wave 
function (6) is written as 



A. Derivation of the amplitude equations 

The stationary condition SS/6Cj = in Eq. (12) re- 
sults in a formal expression of the amplitude equations 



($,|(^|-i/)|vI/)=0. 



(13) 
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(11) 



Here, the time-derivative of the right-wave function is 
decomposed into two parts: 



d_ 
dt 
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(14) 
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with 



Q-spaces, 



i? = ^<clc„ 



'Iq^p'^q: 



(15) 



and rjP = {(f>p\(j)q). This anti-hcrmitian matrix rj^ plays an 
important role in the formulation of the orbital equations 
as we discuss in Sec. III. We insert Eq. (14) into Eq. (13) 
and obtain 



iCi + {<^i\{iD - H)\^) =0. 



(16) 



Furthermore, we substitute Eqs. (3) and (15) into 
Eq. (16) and derive after some algebra the explicit form 

id I = ^sgn(T)C^(^,)(/ij- - itj'j) 

+ lj2^Sn(^)^riifl)^f^^ (17) 



ijkl 



tJ, 



with Cj, = ($/|cJcj|*), Cpi = ($/|cJc[c/Cj|^), and r a 

i ik 

permutation map rearranging strings of orbital indices to 
the ascending order with the sign defined by sgn(r) = 1 
(—1) when T is even (odd). The amplitude equations of 
Eq. (17) are exactly the same as those of the MCTDHF 
theory. In the MCTDHF theory, one can choose all of 
the 77j- to zero and thereby make the amplitude equations 
easier to solve (see, e.g., Ref. [29]). As will be shown in 
Sec. Ill, in the TD-RASSCF theory aU the ryj can not be 
set to zero due to the truncation of the Cl-expansion, and 
one thus needs a more careful simultaneous optimization 
of the expansion coefficients and the orbitals. 



B. Derivation of the orbital equations 

The other stationary condition SS/{S(f>i\ = in 
Eq. (12) using Eq. (14) yields the set of equations to 
be solved for the orbitals 



q \ /sVras 



iCi\<Pi 



s) = 0, 



(18) 



where the one-particle- one-hole state {^1\ = (^|c|cq is 
introduced. The orbital equations need to be solved in 
both the V- and Q-spaces as indicated by the use of the 
index q. Defining projection operators onto the V- and 
Q-spaces by 



2M 



p = Y^\4>,){ct^,\, 



i=l 
1-P, 



(19) 
(20) 



respectively, the time derivative of each orbital is decom- 
posed into two parts, i.e., contributions from the V- and 



(P + QM 

2M 

Ei-^^)^'- 



Q-space orbital equations 



(21) 



One can obtain a formal expression of the Q-space 
orbital equations by multiplying Eq. (18) by a virtual 
orbital bra- vector {(j)a\ from the left, and by using the 
orthogonality between the active and virtual orbitals. 



(■^fHiD - mm = 0. 



(22) 



Equation (22) is a generalization of Brillouin's theorem 
[52] to time-dependent problems. Substituting Eqs. (3) 
and (15) into Eq. (22) and perfoming some algebra with 
the help of Wick's theorem [51], 

cUaclcq = 44>^qCa + S^cj'^q^ (23) 

C^CaCpCj^CgCq — C^CpC^CsCqCa \ O^C^C^CsCq — u^C^C^CsCq., 

(24) 



the Q-space orbital equations read 

E(^^"-^>HE</'t 

with the density matrices p^ and p^. defined by 



(25) 



(vPlcIc,!*) = ^ sgn(r)C;(,pC,, (26) 
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and 



,J' - 



Pzk 



;*|clciQc,|vI') = Y. sgn(r)C;(,..)C/. (27) 



leVuA 



To circumvent explicit numerical treatments of the vir- 
tual orbitals, we use the Q-space projection operator and 
express Eq. (25) as 



^T.Q\'^^^P^ = E^M*)!'^^-)/'' +EQ^'l'^^>'i' 



jki 



(28) 



where the mean-field operator is defined by 

Wtir)^ f 4>l{z')v{r,r')<t>i{z')dz'. (29) 



We now arrive at formally the same Q-space orbital equa- 
tions as in the MCTDHF theory (see, e.g., Eq. (12) in 
Ref. [30]). The density matrices in Eqs. (26) and (27) 
are, however, now calculated based on the RAS scheme. 



2. V-space orbital equations 

We obtain a set of equations for the T'-space orbitals 
when we multiply Eq. (18) by an active orbital bra- vector 
{(j)j I from the left, 



{m{iD-H)m+i 



-feVRAs 



■^l\^i)Ci + e)^Q.m 



{C/I/uVrasj the values need to be known beforehand to 
prepare the values of 77^, , because 77^, is considered to be 

a function of {C'/}/gVRAs ^i^- Pi' ■ One way to easily cir- 
cumvent the use of implicit integrations is by taking only 
into consideration even occupation numbers in the V2- 
space, which removes p^, , and thus the T'-space orbital 
equations result in 



{til, ~ ^Vl' )^k"i' + / A^fci 



ri'm 



M 



tiplier e*. Similarly, from the stationary condition 
6S/\S(f>j) — 0, or by taking the hermitian conjugate of 
Eq. (30) followed by exchanging i and j, we arrive at 
equations containing the same multiplier 



{^\{iD-H)m)-i 



E 

/eVRAs 



c;($/l*;-) + e;. = 0.(31) 



Subtraction of Eq. (30) from Eq. (31) removes the multi- 
plier and gives the T'-space orbital equations or the time- 
dependent Brillouin's theorem for the active orbitals 

in^D - HWj) - (:i>i\{tD - H)\^) = zfi, (32) 

where the time-derivative of the density matrix is 



P^ 



E 

/6Vras 



cn^im 



{n\^i)Ci 



(33) 



In some cases, a set of 77^ 's may be obtained by solving 
Eq. (32). In the MCTDHF theory, however, it is weU 
known that Eq. (32) does not need to be solved, or in- 
deed can not, and the rjl 's are thus often chosen to be zero 
[42, 43]. Such freedom does not exist in the TD-RASSCF 
theory because the Fock space Vras does not consist of 
all possible configurations. It is always possible, however, 
to set 77* = within each subspace Vk {K = 0, 1, and 2), 
i.e., when two orbitals (j)i{t) and (t)j{t) belong to the same 
subspace (see Fig. 1). In the TD-RASSCF theory, gener- 
ally Eq. (32) needs to be solved for combinations («',/') 
to determine rjf, ( = — (77*,, ) j , where indices with sin- 
gle and double prime hereafter mean that the orbitals 
labeled by them belong to different subspaces. Substi- 
tuting Eqs. (3) and (15) into Eq. (32) and computing 
commutators, the explicit expression of the T'-space or- 
bital equations reads 



J2i^4'^^vf)AiC+J2(-h 



m kl 



,kl 



^ipi 



k"l' 



where 



klm 



(34) 
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[Cj, Cj„ , Cf,, 



ci']m 
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To carry out time propagation of the wave packet, the 
set of equations of motion, i.e., Eqs. (17), (28), and (34), 
need to be solved. Notice that the numerical integration 
of the amplitude and the P-space orbital equations re- 
quires elaborate implicit integration schemes: for solving 
the amplitude equations (17) to compute the values of 



k"V 



klm 



pin = o. 

(36) 



The amplitude and the P-space orbital equations are now 
separable and can be easily solved by usual explicit in- 
tegration algorithms. In this work, all the calculations 
are based on Eq. (36), by which all the singly-excited 
configurations are abandoned. However, it should be 
noted that even within this scheme, the wave packet 
partly includes single electron excitation processes due 
to the time-dependent Brillouin's theorem [Eqs. (22) and 
(32)]. In other words, by solving the Q- and T'-space 
orbital equations for a given set of the P-space orbitals 
{0i(t)}f^, we obtain a new set of orbitals {(l)i{t+dt)}f^{, 
which are variationally optimized to take into account the 
single electron processes between the V- and Q-spaces 
and among the Vq-, Vi-, and 7'2-spaces at any instant of 
time t. 



III. 



PARAMETRIC REDUNDANCY 



A. "P-space orbital equations revisited 

In the preceding section, the derivation of the T'-space 
orbital equations was briefly sketched followed by a dis- 
cussion of how to solve them. We now revisit certain 
details of the derivation. Consider the substitution of 
Eq. (33) into Eq. (32): by using a formal expression of 
the amplitude equations (16), we arrive at another ex- 
pression of the T'-space orbital equations 



{^\{tD - H){1 - Il)\m) - (vl/^Kl - Il){tD - i/)!*) 



with a projection operator defined by 

n= ^ !$/)($/!. 

-feVRAs 



= 0, 

(37) 



(38) 



In the MCTDHF theory, i.e., when Vras is replaced by 
Vfci defined above Eq. (2), since the Fock space Vpci 
includes all possible configurations, the left hand side 
of the expression (37) is zero because (1 — n)|'I'*) = 

(^^1(1 — n) = for any combination of i and j. Hence, 
the T'-space orbitals are completely undetermined, and 
one can therefore choose arbitrary anti-hermitian matri- 
ces for 77^ (see, e.g., Refs. [42, 43]). This fact stems from 
the non-uniqueness of {C'/}/eVFci ^^'^ {4'i}'i=i- As is well 



known in the time-independent SCF theory, a unitary 
transformation of the orbitals 



Here it is important to notice that p^, 



when ei- 



E i'^^-)^^-^ 



(39) 



together with the transformation of the expansion coef- 
ficients 

G/ = E---EG.uvGr„U.GS (40) 

keeps the wave function invariant. This property at a cer- 
tain moment of time is called parametric redundancy in 
time- independent quantum chemistry [52]. In the time- 
dependent formulation, however, it is of importance as 
well to consider the time evolution of the unitary trans- 
formation: 

vi = {4>Mi) 

kl kl 

= E GllGk^ + E Gi^'Gk^r,^, (41) 



kl 



which is formally solved in matrix form 



G{t) = .^cxp 



l\rj{t')-rj'{t'))dt' + A 



(42) 



where A is a constant anti-hermitian matrix and £^ de- 
notes time-ordering. This equation means that between 
any two sets of orbitals, there exists a unitary transfor- 
mation at any moment of time and therefore it ensures 
a unique description of the time-dependent dynamics by 
using an arbitrary set of orbitals. This kind of geomet- 
rical structure is an advanced concept of the parametric 
redundancy and is known as principal bundle, in which 
the gauge map defined by ri{t) exists [47, 49]. Exploit- 
ing the gauge degree of freedom, usually the gauge is 
fixed such that 77^ = at all times to make the T'-space 
orbital equations vanish and simplify the system of dif- 
ferential equations. Another useful gauge choice is em- 
ploying rjf = —ihl by which the use of a larger time step 
is sometimes allowed in time propagation [29, 42, 43]. 

In the TD-RASSCF theory, the expression (37) is still 
a trivial identity for combinations {i,j) belonging to the 
same subspace in V, i.e., if (j)i{t) G Vk and (pjit) e Vk 
(^K = 0, 1, and 2). The time-dependent orbital rotations 
within each subspace are hence undetermined. How- 
ever, this is not the case for the combinations {i',j"), 
i.e., if (t)i,{t) G Vk' and (j)j„{t) e Vk" for K' ^ K" 
{K' , K" = 0, 1, and 2). This is because the unitary trans- 
formations (39) and (40) can now be carried out only 
within each subspace, and Eq. (42) is still true in each 
subspace. Therefore, fixing the gauge such that ry^ = is 
satisfied, what remains is to determine the off-diagonal 
block elements 77^, ( = —{rf-n) ) by solving Eq. (34). 



ther i' or j" denotes the index of a T'o-space orbital. Fur- 
thermore, by taking only into consideration even occupa- 
tion numbers in the 7'2-space, p^, and fy'^, vanish, and 
the T'-space orbital equations thereby result in Eq. (36). 
The amplitude and the 7^-space orbital equations are now 
separable. Another prescription to set p'^, —pi, — Q is 
forbiding electron transitions between Vi- and 7'2-spaces. 
This complete-active-space (CAS) scheme [52, 55] gives 
us formally the same T'-space orbital equations [56] . 



B. Related works 

It is informative to briefly discuss two related works. 
The first one is the MCTDH method with selected con- 
figurations (S-MCTDH) [48]. To simphfy the problem, 
the S-MCTDH method ignores the treatment of the V- 
space orbital equations, which are thus assumed to be 
always satisfied, i.e., supposed to be identities, not equa- 
tions. Although the S-MCTDH method works efficiently 
for computing absorption spectra of a pyrazine molecule, 
the method exhibits numerical instability as well for some 
configuration selections conceivably due to the discard of 
the P-space orbital equations. The other related method 
is based on the MCTDHF theory with a truncation of 
the expansion [57]. To reduce the numerical cost, the 
method employs the time-independent expansion coeffi- 
cients, i.e., fixed values throughout the time propagation. 
These two works abandon either the amplitude or the 
7^-space orbital equations as an additional approxima- 
tion, which lowers the accuracy of the description of the 
dynamics. Both the amplitude and the 7^-space orbital 
equations are exactly treated in the present TD-RASSCF 
theory, in which, the P-space orbital equations are simple 
to solve and the computational cost is largely reduced by 
the RAS scheme. Before closing this section, we empha- 
size that the gauge degree of freedom due to the principal 
bundle structure is a key concept behind the treatment 
of the T'-space orbital equations. One can find a discus- 
sion of this issue in terms of the OATDCC theory in the 
supplementary material of Ref . [50] . 



IV. NUMERICAL APPLICATION 

A. Ground state wave function 

We investigate TVe-electron atoms to demonstrate the 
computational efficiency and analyze the convergence 
property of the TD-RASSCF theory by proof-of-principlc 
calculations. The atoms are modeled by ID systems with 
soft-core Coulomb potential: The one-body operator in 
Eq. (4) is 



h{x) 



1^ 
2l^ 



V{x) 



(43) 



Q-space 



"P-space 

{V = Vi® V2) 



Two-electron 
excitation 



±t 



# of spatial-orbitals in V 

7^2-space 

Ma: 

# of spatial-orbitals in V2 



"Pi-space 

Mi: 

# of spatial-orbitals in Vi 



FIG. 2. Illustration of the single-particle Hilbert space used in 
the calculations for the ID beryllium atom. In this case, the 
"P-space is simply decomposed into two active spaces: only 
two-electron transitions between the Vi- and 'P2-spaces are 
now permitted. The numbers of spatial-orbitals in the Vi- 
and 'P2-spaces are expressed by Mi(> 1) and M2(> 0), re- 
spectively, and the total number by M (= Mi -|- M2). In this 
illustrative example of (M, Mi) = (8,3), the electrons are in 
the lowest energy configuration in the "Pi-space. Notice the 
special case Mi = M, where the P2-space disappears and the 
present RAS scheme thereby becomes the MCTDHF theory. 
In the application to the ID helium atom, the same partition- 
ing in the P-space was used. 



with 



V{x) 



V^ 



I 



(44) 



where Z = Nc = 2 for describing a helium atom [58-60] 
and Z ^ Nc — 4 for a beryllium atom [61, 62]. The 
two-body operator in Eq. (5) is 



i;(a;i,X2) 



1 



V(a^l - a^2)2 -f- 1 



(45) 



In this section, the RAS scheme is simplified by elim- 
inating the inactive-core space Vo as shown in Fig. 2: 
only two-electron transitions are allowed between the Pi- 
and Pa-spaces, in which the numbers of spatial-orbitals 
are Mi(> 1) and M2(> 0), respectively, and the total 
number is M(= Mi -|- M2). In this scheme, the TD- 
RASSCF theory is equivalent to the MCTDHF theory 
when Ml = M. 

Table I lists the ground state energy of the ID beryl- 
lium atom for different combinations of (M, Mi). The 
results calculated by the MCTDHF method (Mi = 
M) exactly agree with those given in Ref. [61]. All 
the results are obtained from imaginary time relaxation 
in a box [—25, 25] discretized by the discrete-variable- 
representation (DVR) [63] quadrature points, A'^dvr — 



256, associated with Fourier basis functions. The inte- 
ger in the parenthesis below each energy value gives the 
number of configurations. Figure 3 depicts some selected 
results for the spin-averaged two-electron density 



P2{X1,X2) ^ ^^Pfk d(Ti dcT: 

ijkl 



x||0l(zi)4(z2)||||0,(zi)0K^2)||, (46) 

where [j • • • ]] means the normalized Slater determinant 
(for two-electron systems, p2ixi,X2) is equivalent to the 
absolute square of the spatial wave function). 

In Table I and Fig. 3, for a fixed Mi, the larger M, 
i.e., the more active orbitals, variationally the more ac- 
curate a result is obtained. On the other hand, for a 
fixed M, the use of larger Mi docs not necessarily gives 
more accurate results. To clarify the physics behind this 
convergence behavior, consider how the RAS scheme is 
working: the four electrons are firstly distributed in all 
possible manners in the Pi-space, from which all possi- 
ble two-electron transitions to the Pa-space take place, 
as shown in Fig. 4, where typical configurations realized 
for (M,Mi) = (8,1), (8,2), and (8,6) are illustrated. 
Notice that, relative to the lowest energy configuration, 
the singly-excited configurations are realized only for 
(M, Ml) = (8, 1). The wave function of (M, Mi) = (8, 2), 
however, includes doubly-excited configurations hereafter 
called quasi-singly-excited configurations, in which one of 
the two excited electrons remains in a low-energy orbital 
near the nucleus but the other occupies a high-energy 
orbital far away from the nucleus. These singly- and 
quasi-singly-excited configurations are excluded from the 
method of (M, Mi) = (8,6), which has an advantage to 
take into account collective four-electron correlated con- 
figurations near the nucleus. In short, the larger Mi, i.e., 
the more upward the partition shifts, the more configu- 
rations in the Pi-space but the less configurations in the 
P2-space the wave function includes. 

Look at the line of Mi = 2 in Table I, where the num- 
ber of configurations is largely reduced, and thereby mak- 
ing the computation time shorter. In this line a slow con- 
vergence with respect to M can be seen: starting from 
the Hartree-Fock (HF) energy, the energy value decreases 
and eventually becomes —6.784533 at Mi = 12, which 
still differs from the value —6.785077 obtained from the 
(M, Ml) = (12,12) calculation. The slow convergence 
property is more pronounced in the calculations using 
Ml — 1. This is due to the lack of the collective four- 
electron correlated configurations for describing the com- 
plex electron correlation near the nucleus. On the other 
hand, for a fixed M{> 6), an energy value calculated 
with Ml = 6 is always more accurate than correspond- 
ing ones with Mi = 1 and 2. However, against the su- 
periority of the use of Mi := 6 to Mi = 1 and 2 for 
calculating the ground state energy. Fig. 3 indicates an 
opposite view; the uses of Mi — 1 and 2 are seemingly 
superior to Mi = 6 for a more accurate description of the 
two-electron density in logarithmic scale. In a region far 



TABLE I. Ground state energy in atomic units of the ID beryllium atom {Z — Nc — 4) for different combinations of {M, Mi) 
(see caption of Fig. 2). The integers in parentheses below each energy value show the number of configurations used in the 
calculation. 



M 


2 


3 


4 


6 


8 


10 


12 


Ml = 1 




-6.771296 


-6.775354 


-6.776631 


-6.776764 


-6.776780 


-6.776782 


2 


-6.739450 


(5) 
-6.771296 


(18) 
-6.779805 


(125) 
-6.784224 


(490) 
-6.784501 


(1377) 
-6.784529 


(3146) 
-6.784533 


3 


(1) 


(5) 
-6.771296 


(19) 
-6.775314 


(77) 
-6.779301 


(175) 
-6.779648 


(313) 
-6.779683 


(491) 
-6.779687 


4 




(9) 


(18) 
-6.780026 


(108) 
-6.781293 


(294) 
-6.781591 


(576) 
-6.781627 


(954) 
-6.781633 


6 






(36) 


(112) 
-6.784736 


(364) 
-6.784814 


(792) 
-6.784838 


(1396) 
-6.784843 


8 








(225) 


(399) 
-6.785041 


(981) 
-6.785049 


(1971) 
-6.785050 


10 










(784) 


(1096) 
-6.785072 


(2144) 
-6.785074 


12 












(2025) 


(2515) 

-6.785077 

(4356) 



from the nucleus and especially around xi ~ 2:2, the den- 
sity is remarkably well described in the calculations with 
Ml = 1 and 2. This observation so far interestingly in- 
dicates that taking into consideration the collective four- 
electron correlated configurations is crucial for obtain- 
ing an accurate ground state energy, while the singly- 
or quasi-singly-excited configurations are important for 
the details of the electron density in the region far from 
the nucleus, which is where the tunneling ionization by 
strong lasers takes place, as will be discussed in the next 
subsection. 



Finally we consider an application to the ID helium 
atom {Z = Ng = 2). Under the same numerical condi- 
tion, a direct solution of the Schrodinger equation pro- 
vides the exact ground state energy —2.238257. On the 
other hand, rapid convergence of the ground state energy 
is observed in the MCTDHF calculation with increasing 
M; starting from the HP energy —2.224210, almost con- 
verged energy is already obtained with M — 7 (see also 
Ref. [58] , which lists the ground state energy of the same 
ID helium atom for several values of M). Here it should 
be noted that, for two-electron systems, the wave func- 
tion in the present TD-RASSCF scheme is invariant with 
respect to the position of the partition (see Appendix A). 
This is also the case for general A^'c-electron atoms when 
M — Nc/2+1, because there are two holes which play the 
same role as the two electrons in two-electron systems. 
In the case of the ID beryllium atom with M = 3, the 
TD-RASSCF calculations using Mi = 1, 2, and 3 hence 
provide equal ground state energy as shown in Table I. 



B. High-order harmonic generation 

To investigate the performance of the TD-RASSCF 
theory in a truely time-dependent setting, we consider 
the dynamics of the ID beryllium model atom {Z — Nc — 
4) interacting with a few-cycle near-infrared laser field. 
The effect of the laser is described by adding the dipole 
interaction in the length gauge to the one-body Hamil- 
tonian (43) as 



h{x,t) 



1 (f 



2dx^ 
with the laser field expressed by 



+ V{x)+xF{t)-iW{x), (47) 



F{t) = Fo cos- 



,2 



sinwi, (-r/2 < t<T/2).(48) 



Here Fq is the electric field strength, w the angular fre- 
quency. All the calculations in this section were carried 
out using a larger box [—300, 300(= L)] discretized by 
^DVR = 2048 points. The real time propagation was 
implemented by the fourth-order Runge-Kutta method 
with time step At = 0.005. To cure the electron re- 
flections at the edges of the box, Eq. (47) includes the 
complex absorbing potential (CAP) function defined by 

W{x) = 1 - cos|7r(|x| - Xcap)/[2(L - a;cap)] | with 

s^cap — 250 for |a;| > Xcap and zero otherwise [64]. To 
keep the calculations stable, a further numerical tech- 
nique is needed. In the Q-space orbital equations (28) 
the density matrix is regularized by the following substi- 
tution to prevent it from being singular 



Prcg = p + e cxp ( - p/e) 



(49) 



with a small constants e = 10 ^ [42]. The same regular- 
ization method was used for the tensor AJ,-, in the V- 
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FIG. 3. (Color online) Logarithmic contour plot of the spin-averaged two-electron density p{xi,X2) [Eq. (46)] of the ID 
beryllium atom in the ground state for different combinations of (M, Mi) as indicated on each panel. The leftmost result of 
(M, Ml) = (2,2) corresponds to the HF result. For comparison, all panels include the same dashed (red) representing the 
density calculated with (M, Mi) = (12, 12). Contours differ by a factor of 10, indicating 0.1 for the line of the innermost island. 



space orbital equations (36). The validity of the numer- 
ical data shown below was checked by carrying out the 
same calculations with using larger boxes, denser DVR 
quadrature points, smaller time steps, and different val- 
ues of the CAP parameter, Xcap- 

As an important observable, the HHG spectra are cal- 
culated from the dipole moment in the acceleration form 



s{n) 



N^ 



"^/2 K=l 



dXf, 



(50) 



which is supposed to be more favorable than in the length 
form especially when a CAP function is used [15, 16]. 
There are also other superiorities for the use of the accel- 
eration form to the length form as discussed in Ref. [21]. 



Note that, in Eq. (50), the laser electric field is excluded 
since it does not contribute to the HHG spectrum. 

Figure 5 displays the HHG spectra of the ID beryllium 
atom induced by a laser pulse specified by Fo — 0.0755 
(2.0 X 1014 Wcm-2), oj = 0.0570 (800 nm), and T = 331 
(3 cycles). For comparison, all panels include the same 
dashed (red) line representing the MCTDHF result cal- 
culated with (M, Afi) = (12,12), which is supposed to 
include most electron correlation. In all cases a dou- 
ble plateau structure appears due to the many-electron 
effect and the use of the short pulse. The vertical dot- 
ted lines in Fig. 5 indicate the positions of the first and 
second cutoffs estimated based on the three-step model 
by solving the classical equations of motion for a free 
electron in the laser field as follows: Within the second 
laser cycle, a liberated electron returns back to the parent 
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FIG. 4. Concept of the excited configurations in the TD-RASSCF theory for the ID beryllium atom (see also Fig. 3). These 
are typical excited configurations for (M, Mi) = (8, 1), (8,2), and (8,6) from left to right, respectively. Relative to the lowest 
energy configuration, the singly-excited configurations are realized only for (M, Mi) = (8, 1). The calculation of (M, Mi) — (8, 2) 
includes the doubly-excited configurations called quasi-singly-excited configurations, where, in the 7'2-space, one of the two 
excited electrons remains in a low-energy orbital near the nucleus but the other occupies a high-energy orbital far away from 
the nucleus. These configurations are excluded from the method of (M, Mi) — (8,6), which, however, includes collective 
four-electron correlated configurations near the nucleus. 



ion with the maxinium kinetic energy 3.15C/p, which ac- 
counts for the first cutoff as 3.15t/p + /p — 29. Sw. Here 
the first ionization potential is estimated from the high- 
est occupied orbital energy — 0.3127982(= — /p ) in the 
HF approximation, and Up = Fq/{4uj^) = 0.439 is the 
ponderomotive potential. Slightly after this moment of 
time, another electron already ejected in the first laser cy- 
cle returns to the parent ion. Treating the two electrons 
coherently with neglecting the electron repulsion [65] , the 
nonsequential double recombination emits a photon hav- 



ing the maximum energy 5.04C/p -I- /] 



(1) 



r(2) 



58.7a;, 



which could roughly explain the second cutoff. Here the 
sum of the first and second ionization potentials is esti- 



mated to be 4^^ +4^^ = E'^+-Eg = 1.129997, where E^ 



and i^g"*" are the ground state energies of the ID beryl- 
lium atom and its dication Be^+, respectively, obtained 
by the HF calculations. 

We start the discussion of the spectra in Fig. 5 by inves- 
tigating the structure of the first plateau (0 < fl/oj < 30). 
Although the overall shape of the HHG spectra is sim- 
ilar in all calculations, the TDHF result, i.e., the re- 
suh of (M, Ml) = (2,2) in Fig. 5 (a) shows a signifi- 
cant disagreement with the result of {M, Mi) — (12, 12). 
This failure is because the creation of the first plateau 
is mainly due to the single ionization and recombina- 
tion processes, which are not included explicitly in the 
TDHF wave function. The increase of M removes this 
shortcoming, and the variational improvement of the ac- 
curacy is thus apparent in Figs. 5 (b) and (c), in which 
all the results are in reasonable agreement. Here the 



agreement especially among the (M, Afi) = (6, 1), (6,2), 
(12, 1), (12, 2), and (12, 12) results importantly indicates 
that the singly- and quasi-singly-excited configurations 
play a leading role to reproduce the first plateau. Around 
15 < fl/u! < 30, however, all the results show small dis- 
agreements, which require further analysis beyond the 
scope of the present work. 

Next look into the second plateau (30 < fi/o; < 60). 
The TDHF result in Fig. 5 (a) differs much from the re- 
sult of (M, Ml) = (12, 12). As discussed so far, the non- 
sequential double recombination roughly estimates the 
creation of the second plateau. Although the TDHF wave 
function takes into account the double continua, the two 
liberated electrons are always in the same spatial orbital, 
which results in the poor accuracy of the TDHF method. 
On the other hand, in Fig. 5 (b) and (c), the second 
plateau is roughly reproduced by all methods, includ- 
ing the (M, Ml) = (6,2) and (12,2) methods despite of 
the large reduction of the number of configurations in 
these approach. It seems to be a formidable task to ex- 
actly reproduce the fine structure. The rich structure in 
the second plateau especially around the second cutoff 
is due to the interference among several quantum tra- 
jectories of the two electrons coming back to the parent 
ion [65]. It is thus still questionable whether the con- 
vergence is completely achieved even in the calculation 
with (M,Mi) = (12, 12). Notice that, in Ref. [14], HHG 
spectra were calculated for a four-electron molecule us- 
ing a similar laser pulse by the MCTDHF method. In 
this related work, however, all the results including the 
TDHF one show reasonable agreement in both the first 
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FIG. 5. (Color online) HHG spectra of the ID beryllium atom calculated for different combinations of (M, Mi) (see caption of 
Fig. 2). The result of (M, Mi) = (2,2) shown by the solid black line in (a) corresponds to the TDHF result. For comparison, 
all panels include the same dashed (red) line corresponding to the result of (M, Mi) — (12, 12). The laser pulse is specified by 
the parameters: Fo = 0.0755 (2.0 x lO" Wcm"^), uj = 0.0570 (800 nm), and T = 331 (3 cycles). The first and second cutoff 
energies are estimated to be 29.8aj and 58.7aj, respectively, as shown by vertical dotted lines (see text). 



and second plateaus. The convergence behavior of the 
MCTDHF calculations will thus sensitively depend on 
the system, as well as the parameters of the driving laser. 



V. CONCLUSIONS AND OUTLOOKS 

In this paper, a new theoretical framework for de- 
scribing time-dependent many-electron dynamics, de- 
noted TD-RASSCF theory, was proposed and formulated 
on the basis of the time-dependent variational princi- 
ple. The key concepts of the theory are the use of the 
time-dependent orbitals and the truncation of the CI- 
expansion by incorporating the RAS scheme. Abandon- 
ing the full-CI expansion gives rise to important changes 



in the formulation as compared to the MCTDHF theory. 
The TD-RASSCF theory thereby requires us to solve 
the T'-space orbital equations. To make the amplitude 
and the P-space orbital equations separable, we only al- 
low transitions of even numbers of electrons between the 
Vi- and ■P2-spaces. In a proof-of-principle application 
to the ID beryllium atom, the TD-RASSCF method ex- 
hibited a reasonable convergence behavior with accumu- 
lating the number of active orbitals in both calculations 
of the ground state wave function and the HHG spectra 
induced by an intense laser pulse. By shifting the posi- 
tion of the partition between the two active spaces and 
changing the number of active orbitals in each subspace, 
one can flexibly take into account the electron correlation 
needed for describing the phenomena of interest. This 



flexibility and the accompanying gain in computational 
efforts allow us to promote the TD-RASSCF theory as 
a promising method for application in larger atoms and 
molecules, beyond the reach of methods based on full CI 
expansions. 
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means the rest of the Slater de- 
Defining two new orbitals by 

*M-ita)) - cos0\(j>M~it(i)) - sin6'|(/>Mta)> ^^d 
= sin6'|0M-it(i)> + cos6l|(/)Mt(i)> such that 



where |$) 
terminants. 



l'^Mt(i) 

the condition C'm-i sin 6* + Cm cos 9 = is fulfilled 
instance, the wave function is rewritten as 



for 
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Appendix A: Wave function of two-electron systems 

Consider a two-electron system in the MCTDHF 
method with M(> 2) active spatial-orbitals. Ex- 
pressing the index of the active orbitals \(f>i) by i — 
It, li, • • • , Af ti MX, the wave function reads 



I*) = \^)+CM-i\\(pit9M-n} - \(pn9M-it} 

+Cm { \4>i-t4'Mi) - |'/'u0Aft) ) 1 (Ai) 



I*) 



m + Cm 



M-1\\V11-V'M-1I/ 



'^nvM-if/ h' 



,(A2) 



where C^,f_i = Cm-^i cos 6* — Cm sin0. Proceeding with 
the orbital manipulations, one can eliminate any or even 
all of singly-excited configurations relative to the lowest 
energy configuration. In the TD-RASSCF scheme, all the 
singly-excited configurations from Vi to V2 are likewise 
removable. Thus the wave function in the TD-RASSCF 
theory is invariant with respect to the position of the par- 
tition between Vi and V2- Finally note that, exploiting 
this flexible structure of the two-electron wave function, 
one can ultimately arrive at the expression of the wave 
function in terms of geminals [66] instead of orbitals. 
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